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ABSTRACT :  The  approximation  of  a  function  by  a  sum  of  complex 
exponentials,  in  general,  is  a  nonlinear  optimization  problem. 
The  optimization  problem,  however,  is  linearized  through  the 
application  of  the  pencil  of  function  method.  This  non¬ 
iterative  method  yields  the  best  exponential  approximation 
for  a  given  order  of  approximation.  The  method  differs 
radically  from  the  classical  Wiener  least  squares  approach 
in  the  sense  that  exponents  calculated  by  the  pencil  of 
function  method  are  directly  proportional  to  the  integrated 
squared  error  in  the  approximation.  As  the  integrated  squared 
error  approaches  zero,  the  exponents  calculated  by  the  pencil 
of  function  method  approach  the  best  least  squares  exponents 
in  a  continuous  fashion.  Among  the  advantages  of  the  method 
are  its  natural  insensitivity  to  noise  in  the  data  and 
explicit  determination  of  the  signal  order.  Examples  are 
presented  to  illustrate  the  stability  of  this  technique 
especially  when  noise  is  present  in  the  data. 

»\ 


This  work  has  been  supported  by  the  Office  of  Naval  Research  under 
contract  N00014-79-C-0598. 


1.  Dept,  of  Electrical  Eng.  2. 

Rochester  Inst,  of  Tech. 

Rochester,  NY  14623 

3.  Department  of  Electrical  Engineering 
University  of  South  Florida 
Tampa,  Florida  33620 


Dept,  of  Electrical  Eng 
Syracuse  University 
Syracuse,  NY  13210 


Distribute  oo  Unbmited  __ 

81  11  12  077 


UNCLASSIFIED 


AECUftlTY  CLASSIFICATION  OF  THIS  PACE  (Whmn  Dm  la  Entarod) 


REPORT  DOCUMENTATION  PAGE 


.  REFORT  NUMBER 


TR-81-8 


4.  title  (and  Submit) 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


S.  RECIPIENT'S  CATALOS  NUMBER 


s.  type  of  report  •  period  covered 


Approximation  by  a  Sum  of  Complex  Exponentials  Technical  Reoort  No.  8 
Utilizing  the  Pencil  of  Function  Method 


S.  PERFORMING  ORG.  REPORT  NUMBER 


7.  AUTHORf*; 


T.K.Sarkar,  D.D. Weiner  and  V.K.Jain 


B.  CONTRACT  OR  GRANT  NUMBERfaj 


N00014-79-C-0598 


t.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Dept,  of  Electrical  Engineering 
Rochester  Institute  of  Technology 
Rochester,  New  York  14623 


(<■  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Department  of  the  Navy 
Office  of  Naval  Research 
Arlington,  Virginia  22217 


.  MONITORING  AGENCY  NAME  A  ADDRESSffl  dlllaranl  I 


10.  PROGRAM  ELEMENT.  PROJECT,  TASK 
AREA  B  WORK  UNIT  NUMBERS 

61153N  RR021-01 

RR021-01-01 

NR371-014 


IZ.  REPORT  DATE 

November  1981 


IS.  NUMBER  OF  PAGES 

20 


i  Controlling  Oltlca)  IS.  SECURITY  CLASS,  (ol  thla  ration) 

UNCLASSIFIED 


IS.  DISTRIBUTION  STATEMENT  (of  thla  Raport)  _  _ 

— ffisfRjr.u  a 


17.  DISTRIBUTION  STATEMENT  (ol  tha  abatract  antarad  In  Block  SO,  II  dlllarant  bom  Rapori) 


19.  KEY  WOROS  (Contlnua  on  ra  vat  at  a/da  tf  nacaaaary  and  Idantlfy  by  bloc a  mmbar) 


1)  Approximation  Problem 

2)  Complex  Exponentials 

3)  Pencil  of  Function 

4)  Pole-zero  Modeling 

5)  Wiener  least  squares  _  _ 


20.  ABSTRACT  (  Continue  on  rovorao  oldo  It  nocoaaary  a n«f  Identity  by  block  number) 

The  approximation  of  a  function  by  a  sum  of  complex  exponentials,  in  general, 
is  a  nonlinear  optimization  problem.  The  optimization  problem,  however,  is 
linearized  through  the  application  of  the  pencil  of  function  method.  This 
non-iterative  method  yields  the  best  exponential  approximation  for  a  given 
order  of  approximation.  The  method  differs  radically  from  the  classical 
Wiener  least  squares  approach  In  the  sense  that  exponents  calculated  by  the 
pencil  of  function  method  are  directly  proportional  to  the  integrated  squared 
error  in  the  approximation.  As  the  integrated  squared  error  approaches  zer 


form  li7- 

I  JAN  79  W * 


EDITION  OF  I  NOV  aa  is  OBSOLETE 
S/M  01 02- LF -01 4-4401 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  THIS  RAGE  (Whan  Data  I 


SECURITY  CLASSIFICATION  OF  THIS  RACE  (Wi«i  Dim  Mmirntmt) 


20. 

the  exponents  calculated  by  the  pencil  of  function  method  approach  the 
best  least  squares  exponents  in  a  continuous  fashion.  Among  the  advantages 
of  the  method  are  its  natural  Insensitivity  to  noise  in  the  data  and 
explicit  determination  of  the  signal  order.  Examples  are  presented  to 
illustrate  the  stability  of  this  technique  especially  when  noise  is  present 
in  the  data. 


unclassified 


I.  INTRODUCTION 


The  approximation  of  an  arbitrary  real  function  x(t)  which  is 
nonzero  for  t  _>  0,  by  a  sum  of  complex  exponentials  often  leads  to 
a  nonlinear  least  squares  problem.  Recently  Goulb  [1]  has  presented 
a  class  of  least  squares  problems  in  which  the  variables  separate. 
One  of  the  objectives  of  this  paper  is  to  show  that  the  pencil  of 
function  method  also  decouples  the  approximation  problem  [2-3). 

In  this  presentation  we  apply  the  pencil  of  function  method  to  both 
the  continuous  function  x(t)  and  y(t),  which  is  the  discretized 
version  of  x(t).  We  show  that  the  resulting  expressions  are  quite 
different  for  the  two  cases.  First  we  present  the  discrete 
approximation  problem  followed  by  the  continuous  case.  Finally 
examples  are  presented  to  illustrate  the  stability  of  this  method 
in  the  calculation  of  the  exponentials  particularly  for  noisy  data. 

2.  Definition  of  the  Problem 

The  problem  of  interest  is  to  approximate  an  arbitrary  square 
integrable  function  x(t)  by  a  sum  of  complex  exponentials,  i.e. 

£  “l 

x(t)  ~  l  l  A,,  •  (t)k  1  •  exp{s  t}  (1) 

i=l  k=l  1 

where 

l 

l  m  =n  (2) 

1*1  1 

Here  are  referred  to  as  the  poles  of  the  system  and  m^  is  the 
order  of  the  pole  s^  n  is  the  total  number  of  poles  counting 
the  multiplicities  of  repeated  poles. 


(1) 


In  most  problems,  however,  we  do  not  deal  with  continuous 


functions  but  with  sampled  sequences.  So  we  will  assume  that  the 
sequence  (y (p ) }  has  been  sampled  at  a  sufficiently  fast  rate  from 
the  continuous  function  x(t)  so  that  no  information  has  been  lost. 
Thus  for  a  sampled  data  system  we  have 

a  mi 

y(p)  =  l  l  A  .  {p}k_1  (z  }P  (3) 

i=l  k=l  1 

where 

=  exp(s^T)  (4) 

where  T  is  the  sampling  interval.  Here  z^  are  referred  to  as  the 
z-domain  poles. 

The  problem  of  approximating  A  and  s  in  (1)  or  A  anA  z  in 

IK  X  liC  X 

(3)  is  a  nonlinear  approximation  problem.  However  we  can  linearize 
this  nonlinear  problem  through  the  introduction  of  the  T  operator. 

3.  Properties  of  the  T  Operator. 

The  T  operator  Is  defined  to  be  an  operator  similar  to  a  reverse 
integral  operator.  It  is  defined  on  the  space  of  squared  integrable 
functions  in  the  following  way: 

CO 

T[Y]  A  T{y(p)}  A  {  l  y(k)}  (5) 

k=p 

so  for  example 

T{0, 1,0,0, . . . }  =  {1,1, 0,0,0,...}  (6) 

and 

T{p(0.5)P}  =  {2p(0.5)P  +  2(0. 5)P}  (7) 

From  the  definition  of  T  in  (5)  we  can  prove  the  following: 

(a)  T  is  a  linear  operator,  since 

T[aY+gZ]  =  aT[Y]+0T[Z]  (8) 


(2) 


where  a  and  £S  are  scalar  quantities.  [Y]  and  [ Z j  are  two  sampled 
data  sequences  as  defined  in  (5)  . 

(b)  A  constant  multiplication  of  a  unit  impulse  remains  unchanged 
by  the  application  of  the  T  operator,  i.e. 

T[a6]  A  T{a,0,0,...}  «=  aT{l,0,0, . . . }  =  aT[6]  (9) 

(c)  for  simple  poles  (i.e.  m^=l  for  all  i)  we  have 

n  n  (z.)P 

T[Y]  =  T(y(p)}  =  T{  l  A.(z  )P}  =  {  J  A  )  (10) 

i=l  11  i=l  i 


(d)  Successive  operations  of  T  are  obtained  in  the  following  way 


and  thus  we  can  generate  a  polynomial  in  T.  One  such  polynomial. 


(3) 


Furthermore 


[l-(l-z)T]  {y(p)}  *=  (0}  for  k  ^m 


(16) 


Geometrically,  the  above  property  asserts  that  the  operator 

k 

[l-(l-z)T]  projects  Y  into  a  lower  dimensional  hyperplane  in  z of 
order  m-k.  Hence,  it  follows  that,  if  the  operator 


n  [l-(l-z)T]  i 
i=l 

is  applied  to  [Y]  which  is  given  by 

vk-1 


m. 

n  i 


[Y]  =  {y(p)>  =  {  l  l  A.,  (p)K_X(z  )P> 
i=l  k=l  1 


(17) 


(18) 


then  the  result  would  be  the  null  sequence  {0}.  We  have  thus 
linearized  the  problem  by  first  attempting  to  estimate  the  poles  z^. 
After  the  poles  have  been  computed,  the  residues  A^  at  the  poles 
can  be  obtained  by  a  least  squares  procedure. 

5.  Determination  of  the  Poles  z^. 

In  order  to  determine  the  poles  z^,  we  define  the  following 


operations.  Let 


[y2]  A  t[yl]  A  T[TYq]  =  T  [Yq] 


(19) 


and 


[Z±]  A  [l-(l-z)T][Yi]  =  [1+dT] [Y±] 


(20) 


where 


d  A  -(1-z)  (21) 

From  (20)  it  is  seen  that 

[Zi]  =  [Yi]-W[T]tY1]-[T]{Yi_1+d[T][Yi_1]}  *=  [T]  [Z^]  (22) 

The  sequence  [Z^]  is  called  a  pencil  of  functions  parameterized 
by  the  scalar  parameter  d  [2-5],  It  is  a  linear  combination  of  functions 
[Y^]  and  [Y^+^]  through  the  parameter  d.  The  pencil  of  functions 
contain  very  important  characteristics  from  a  system  identification 


W 


r 


point  of  view.  When  the  set  [Z^]  of  n+l  functions  becomes  linearly 
dependent  then  (d+1)  becomes  a  system  pole.  Nov  we  define  the  gram 
matrix  [H]  whose  ith  row  and  jth  column  are  defined  as 
[Hi  -  (h  A  <Zl.Z.>  -  <Il«THY1],U+dTIIT  ]>) 

■  ld*<d<1'i+i-'IJ+i>+<''i-''j+i>w<Yi+rYJ>+<1'1'Yj>»  <23> 

where 

OO 

<y.,y  >  =  l  {y.(p)Hy,(P)}  (24) 

3  p=0  3 

and  d*  denotes  the  complex  conjugate  of  d.  The  underbar  distinguishes 
a  matrix  from  a  sequence. 

The  significance  of  the  gram  matrix  [H]  lies  in  the  fact  that  the 
set  [Z11 ,  [Z  ],  [Z3]  ...  [Zk]  are  linearly  independent  if  and  only  if 
det[H]>0,  whereas  it  is  linearly  dependent  if  and  only  if  det[H]=0. 
Thus  in  order  that  there  are  n  poles  in  the  sequence  [Z^]  it  is 
necessary  and  sufficient  that  det[H]  be  positive  or  zero  according 
as  k<n,  or  k>n.  [This  is  clear  from  (16)] 

After  some  very  tedious  algebraic  manipulations  it  can  be  shown 


that 


Ic.  k. 

det[H]  =  l  l  (d*)1"1(d)j_1Q 
i=l  j-1 


ij 


(25) 


where  are  the  diagonal  cofactors  of  the  matrix  [G]  whose  elements 


are  defined  as 


A  <Yi»Yj>  =  Io{yi(p)}{yj(p)}] 


(26) 


An  on  line  procedure  for  determining  (26)  is  given  by  the  following 
equations : 


xi(k>  *  y±(k) 


Xj(k) 


m-1 

I  x^v) 
v-k 


I 

i 

s 


5 


5ij 


Mi  1 

=  (-l)i+^  l  [x  (k)x  (k)-  l  (k-m)1  Vx  (M)x  (k)-  *  (k-M)^  \ ,,(M)x  (k) 
k=0  2  v-1  v  2  y=l  M  1 

+  I  1  x  (M)xi(M)(k-M)i+j_V“y] 


V=1  V=1 


v  y 


and 


M-l 

xi (M)  =  l  x1_1(v)- 


v=0 


Here  it  is  assumed  that  x^(p)=0  for  p>M.  This  assumption  is  reasonable 
for  practical  problems. 

When  the  order  of  the  approximation  k  in  (25)  becomes  equal  to  the 
order  of  the  system  the  sets  [Z^],  [Z^],  ...  [Z^]  becomes  linearly 
dependent  and  hence 


k  k 

det [H]  =  0  =  l  l  (d*)i_1(d)j_1Q  . 

i=l  j=l  2 


Sines 


Old  +  •  V>JJ 


and 


det[G]  =  0  when  det[H]  =  0, 

(27)  can  be  rewritten  as 

n+1  .  . ,  n+1  .  . , 

(  SJQ.,W*)"-i+1)-(  Z^Q,.(d)“-J+1}  -  0 

i=l  j=l 


(27) 


(28) 


(29) 


(30) 


From  (30)  the  poles  can  be  obtained  from  the  solution  of  the  following 


polynomial  equation: 

n+1  ... 

!JQii(-l+z)n=j+1  -  0 


i=l 


(31) 


and  the  poles  s^  are  obtained  from  through  the  transformation 

-  y  ln  fzj 

Since  the  function  to  be  approximated  is  considered  to  be  real 

* 

there  will  also  be  a  complex  conjugate  set  of  poles  s^  as  indicated  by  (30), 


(6) 


Once  the  poles  or  z^  are  determined  the  unknown  constants  A 
can  be  determined  from  a  least  squares  procedure. 

So  far  we  have  done  the  analysis  for  discrete  systems.  However 
the  analysis  for  the  continuous  system  will  be  done  in  section  6. 

5.  Error  Analysis  in  the  Presence  of  Noise. 

When  the  signal  x(t)  or  the  sequence  (y(p)}  does  not  contain  any 
noise  then  the  poles  z^  given  by  (31)  are  the  exact  poles  of  the 
system.  However  if  noise  is  present  then  they  no  longer  represent  the 
system  poles.  The  objective  of  this  present  section  is  to  find  the 
error  in  the  location  of  the  poles  due  to  additive  noise. 

Let  (e(p)}  be  the  noisy  sequence  and  let  there  exist  an  approximant 
of  (e(p)}  and  call  it  (y(p)},  which  provides  the  best  exponential 
approximation  of  the  signal  in  the  given  subspace.  The  assumption  made 
here  is  that  the  optimum  approximant  is  unique. 

For  convenience  let  us  call 

e(p)  -  y (p)  +  er(p)  (32) 

We  now  introduce  the  error  sequence  in  the  normalized  form 

[ (e(p)  -  y(p) }]  =  e(r(p)}  (33) 

where  e  is  a  non-negative  real  number  chosen  to  establish  the  equality 

£{r(p)}2  =  £{y(p)}2  (34) 

P  P 

Clearly  e  represents  the  fractional  error  in  the  representation  and  is 
therefore  a  measure  of  the  degree  of  approximation  of  (q(p)}.  If 
(y(p)}  is  the  best  approximant  of  (e(p)},  then  (r(p)}  must  be 
orthogonal  to  (y(p)l.  Thus,  in  particular,  e-0  characterizes  the 
perfect  representation. 

Now  consider  the  following  sequences: 


(7) 


(35) 


[E^]  =  {e1(p)}  =  (e(p) } 
[E2]  =  {e2(p)}  =  T{ei(p)} 


[Ek+l3  =  {ek+l(p)}  =  T{ek(p)1  =  T{yk(p)+erk(p)} 

In  (35)  {y (p ) }  and  {r(p)}  are  not  knovm  explicitly , nor  is  e  for  that 
matter.  We  begin  with  the  computation  of  the  (n+l)X(n+l)  gram  matrix 
[B]  for  the  sequence  [E^],  [E2],  [E^],  i-e* 

[B]  =  [b±  A  <{e±(p) } ,  {ej(p)}>] 


”y1(p)  +  er1(p) 


•  Jy1(p)+Er1(p) . yn+l(p)+£rn+l(p) 


Jn+l(p)+ern+l(p)J 

Since  { y  (p ) }  is  a  sequence  of  order  n,  the  corresponding 
[Y^] ,  ....  tyn+1l  sequences  are  linearly  dependent.  Thus  there  exists 
then  a  nonsingular  matrix  [P]  such  that 


[Y]  A 


{y1(p)} 


{yn+l(p)}' 


-  tp  ]  J 


V 


-  [P] [Y] 


(37) 


From  the  theory  of  orthoganal  transformation  the  sequences 
(y1(p)},  (y2(p)},  ...»  (yn(p)}  must  form  an  orthonormal  set. 
Similarly  we  obtain 


tR] 


{rL(p)} 


{rn+l(p)}’ 


(ri(p) 


} 


=  [P] 


{rn+l(p)J 


}  =  [P][R] 


(38) 


Observe  in  (38)  the  sequences  (r^(p)),  ....  {rn+^(p)}  are  not  orthonormal 
to  each  other.  Substitution  of  (37)  and  (38)  into  (36)  yields 


(36) 


(8) 


(39) 


[B]  =  [P][Y+eR][Y+cR]TtP]T 

=  [P] ( [Y] [Y]T+e(R] [Y]T+c [Y] [R]T+c2 [R] [R]T] [P]T 
Observe  that 


[Y][Y)T 


.  ;0, 

(^3  I  O  J(n+l)x(n+l) 


(40) 


where  I  is  the  nxn  identity  matrix  and 

(r1(p)}*{y1(p)}  -  (r1(p)}’{yn(p)} 

{7n+i(p)}-{yi(p)}  ...  {rn+1(p)}-{7n(P)}. 


[R][Y]T  = 


(41) 


By  utilizing  the  theorem  on  determinant  expansion  as  presented  in 
Appendix  A,  we  get 

i 

I  ^  l  e{t  (p)}*{y  (p)} 
nxn  ,  n+1  1 


det [B]  =  (det [P] )  *det 


A  n+r  _ _ _ 

O  'l  e2{7n+l(p))'{yn+l(p)} 

I  * 


+  terms  of  the  order  of  e  and  higher 

£2'M{rn+l(p)1!|2,(detl-1)2  +  0(c3)  + 


(42) 


Thus  it  is  clear  that  yj det[B]  is  of  the  same  order  of  magnitude 
as  the  error  measure  £,  i.e. 


det [B]  =  0(e)  (43) 

In  particular,  a  neccessary  and  sufficient  condition  for  perfect 
approximation  of  a  signal  of  order  n  is  that  det[B]  must  vanish. 

In  order  to  find  the  error  in  the  location  of  the  poles  first  the 
error  in  the  diagonal  cofactors  of  the  gram  matrix  [B]  has  to  be 
estimated.  The  diagonal  cofactors  are  given  by: 


(9) 


Dii  -  det 


y1(p)+£r1(p) } 


J{y1_1(p)+er  (p)}[  j 


V_<WP)+CWP),_/ 


V 


(y1+1(p)«tl+1<p)(-.  •  (yn+1<P)«  Vl(p> 1 


(44) 


Let  [Sj  be  the  matrix  which  orthonormalizes  the  vectors,  from 

{yL(p) }. . . (y^Cp) Hyi(p) } •  •  •  (yn+1(p)  1  to  (y^p) } . . . {y^^p) } (yi+1(p) } 

. . . (y  ,,(p)}.  So  we  have 
n+1 


{y-^pH  \ 


{yi-i(p)} 

{y1+i(P)} 


}  =lV 


V{yn+l(p)}y 


f {yx(p)}  ^ 

<  (y1+1(p>)  >*  'W 


V{y«+i(p)}y 


(45) 


{^(p))  \ 


{ri-l(p)}  V 

tr1+i(p)}  l 


{rn+l(p)} 


[S,] 


f{Tl(p)}  ) 


J(r  (p))l 


V{rn+l<p))y 


(46) 


Observe  that  the  vectors  {r,(p)}  ...  {r  , (p) )  are  not  a  set  of 

1  n+1 

orthonormal  vectors.  Substitution  of  (45)  and  (46)  into  (44)  yields 
D,,  -  det 


'ii 


I' 


(det [S. ]) 


•det 


[Y,  1  [1a  ]T+£[R,  ]  [±,  ]T+e[YinR1IT+e2  [Rj  [R^1]  (47) 


(10) 


By  utilizing 

[y.][y.]t=  [I] 

— i  l  nxn 

and  applying  the  theorem  on  determinant  expansion  [as  presented  in 
Appendix  A]  we  obtain 


n+1 


ii 


Qii  +  2e(det[S1D 


'*<  l  l  {y,(p)>*{r  (p)}1 
1  J-1  P-0  3  J  I 

i^j 


(48) 


where  Q  are  the  diagonal  cofactors  of  the  gram  matrix  formed  by 

the  sets  [Y, ] ,  [Y„]  ...  [Y  , ]  and  is  given  by  (26).  The  underlying 
1  z  n+l 

assumption  in  (48)  is  of  course  we  have  the  right  order  of  the  system 
n,  i.e. 


det [G]  =  0.  (49) 

Observe  that  (48)  implies  continuous  differentiability  of  the  diagonal 

cofactors  with  respect  to  the  error  e,  so  that 

lim  D  =Q  for  all  i.  (50) 

£~*0  11  1 

Let  z^  be  the  approximate  poles  that  are  obtained  by  utilizing  the 
diagonal  cofactors  instead  of  in  (31).  So  the  poles  z^  are 
the  roots  of  the  following  equation 

i/D..(-l+z!)n+1~i  =  0  (51) 

i=l  ^ 

Substitution  of  (48)  into  (51)  illustrates  that  the  approximate  poles 
z^  are  related  to  the  exact  poles  z^  by  the  relationship 

z[  =  Zi  +  +  0(e2)  (52) 

where  are  certain  constants.  The  approximate  poles  z|  are  again  a 
continuous  function  of  6.  So  the  error  in  the  approximation  (z^-z^) 
is  a  continuous  function  of  e  and  approaches  zero  uniformly  as  e-*0. 


(ID 


This  is  the  unique  feature  of  this  technique.  The  continuous 
dependence  of  the  error  in  the  pole  locations  on  e  is  not  guaranteed 
by  other  optimization  techniques  and  the  classical  nonlinear  methods 
based  on  Wiener  filter  theory  [4-5]. 

6.  Continuous  Version  of  the  Pencil  of  Function  Method. 

In  the  previous  sections  the  approximation  problem  involved 
discretizing  a  continuous  function  which  was  then  processed  digitally. 
A  strictly  continuous  version  is  also  available.  The  problem  is  to 
approximate  the  continuous  function  x(t)  in  the  form  given  by  (1),  i.e. 

Z  m± 

x(t)  =  I  I  A  {t}k  1  exp{  s  t} 
i=l  k=l  1 

In  this  case  we  want  to  estimate  s^  instead  of  expts^T}.  For  the 
continuous  case  the  T  operator  of  the  previous  sections  is  replaced 
by  the  reverse  time  integral  operator  so  that 

OO 

Xi+l(t)  =  /xi(f>df  (53) 

and  x^(t)  =  x(t) 

The  two  sided  Laplace  Transform  of  (53)  is  given  by 


i+1 


00  X  (s) 

(S)  =  /  Xi+l(t)e  dt  =  i~ 


(54) 


and 


<•>  -  l  l  V 

i=l  j=l  s  (s-si);] 


N(s) 


n  (s-s  )mi 
i-1  1 


(55) 


Next  we  consider  the  pencil  of  functions  defined  by 

{x1(t)-Ax2(t)}>  (x2(t)-Xx3(t)},  ....  (xn(t)-Axn+1(t)}  (56) 

Since  for  this  problem  it  is  simpler  to  deal  with  the  Laplace 
transform,  attention  is  focused  on  the  set  of  pencil  of  functions  defined  by 
(X1(s)-AX2(s)},  {X2(s)-AX;j(s)} .  {V8>"XXn+l(8)}  (57) 


(12) 


The  above  set  of  pencil  of  functions  displays  some  interesting 
properties  if  checked  for  linear  independence.  Checking  for  linear 
independence,  we  write 

l  [ak{Xk(s)-AXk+i(s)}]  =  0  (58) 

k=l 


Substitution  of  (54)  into  (58)  yields 

|  sn_k_1[a  (s-A)X(s) ]  =  0 
k=l  K 


V  n-k 

l  s 

k=l 


[ak(s-A) 


Z 

n 

i=l 


N(s)  j 
(s-Si)mi 


0 


(59) 


Equation  (59)  represents  a  polynomial  in  s  of  the  order  n.  Since 

1,  s,  s^,  ....  sn  1  span  a  n  dimensional  subspace,  this  set  is  linearly 

independent  and  hence  the  set  of  n  coefficients  a^  must  be  identically 

zero  if  (59)  is  to  be  satisfied  for  all  values  of  s.  The  above  holds 

as  long  as  A  ^  A  .  For  A  =  A^,  it  turns  out  that  (59)  has  a  root 

2  n-2 

s  =  A^  which  can  be  factored  out.  The  set  1,  s,  s  ,  s  spans 

a  (n-1)  dimensional  subspace.  Nevertheless  there  remains  n  set  of 
coefficients  a^.  It  is  now  possible  for  (59)  to  be  satisfied  with 
at  least  one  non-zero  coefficient.  This  implies  linear  dependence  of 
the  set  (57)  when  A  =  A  .  So  if  A  is  a  system  pole  then  the  gram 
determinant  of  the  set  of  functions  in  (57)  will  be  zero.  The  gram 
determinant  of  the  given  set  in  (56)  is  defined  by 

r<Xl-Ax2,  Xl-Xx2>  .....  <Xl~Ax2>  xn“Axn+1>  *"j 


[F]  - 


(60) 


L.<Xn"Axn+l’  XrXx2>  * 


<x  -Ax  , , ,  x  -Ax 
n  n+1  n  n+1 


(13) 


where 

00 

<x  ,  x  >  =  /x  (t)x  (t)dt  (61) 

1  j  o  1  •* 

Thus  det[F]  =  0  if  A  Is  a  system  pole.  After  some  tedious  algebraic 
manipulations  it  can  be  shown  that 


n+1  n+1  ,  ,, 

det[F]  =  ll  v/A  An  H  l  -^/A*  A  ^  )  =  0 

-  -  11  j=l  -1-1 


(62) 


i=l 


where  A^  is  the  kth  diagonal  cofactor  in  the  grammian  [F1]  which  is 
defined  by 


[F1] 


<X1’X1> 


L<Xn+l’  V 


<Xl’Xn+l> 


<Xn+l‘  Vl"  J 


(63) 


The  polynomial  equation  whose  roots  are  the  pole  locations,  is  then 
given  by 

i-i+1 


n+1  i— 


=  0 


(64) 


After  the  poles  have  been  obtained  the  residues  are  obtained  by  a 
least  squares  fashion. 

7.  Examples 

As  an  example  consider  the  transient  response  of  a  conducting  pipe 
tested  at  the  ATHAMAS-I  Electromagnetic  Pulse  (EMP)  simulator.  The 
conducting  pipe  is  10m  long  and  lm  in  diameter.  Hence  the  true  resonance 
of  the  pipe  is  expected  to  be  in  the  neighborhood  of  14MHz.  Also  the 
pipe  has  been  excited  in  such  a  way  that  it  is  reasonable  to  expect  only 
odd  harmonics  at  the  scattered  fields.  The  data  which  have  been 
measured  are  the  integral  of  the  electric  field  and  hence  is  available 
in  terms  of  a  voltage.  Thus,  in  addition  to  the  frequencies  of  the 
conducting  pipe  one  should  also  observe  a  very  dominant  low  frequency 
pole.  The  same  transient  data  as  depicted  in  Figure  1  is  used  for 


(14) 


analysis.  The  results  for  a  fifth  and  a  seventh  order  system  are  as 
follows : 


For  n  *  5,  the  poles  in  radians/sec  are 
(-0.0029  +  JO. 083)  x  109 
(-0.0428  +  JO. 217)  x  109 
(-0.0098  )  x  109 

For  n  =  7,  the  poles  in  radians /sec  are 
(-0.0058  +  j0.084)  x  109 
(-0.0270  +  .10.219)  x  109 
(-0.0270  +  jO.550)  x  109 
(-0.0012  )  x  109 


(=13. 33MHz) 
(=35.20  MHz) 
(=  1.56  MHz) 


(=13.40  MHz) 
(=35.10  MHz) 
(=87.60  MHz) 
(=  0.19  MHz) 


It  is  Interesting  to  observe  that  the  real  pole  due  to  the 
integrator  has  been  obtained.  This  pole  is  a  very  dominant  pole  as 
the  data  were  recorded  after  having  passed  through  an  integrator. 
The  above  results  display  a  dynamic  range  of  approximately  1000:1 
for  the  values  of  poles  of  the  conducting  pipe. 

Next  the  data  were  differentiated  to  get  rid  of  the  undesirable 
dominang  pole  of  the  integrator.  The  differentiation  was  done 
numerically.  For  a  fourth  and  sixth  order  system  the  above  results 
have  been  recalculated  as  follows: 

For  n  =  4,  the  poles  in  radians/sec  are 

(-0.0026  +  j0.086)  x  109  (=13.70  MHz) 

(-0.0480  +  jO. 235)  x  109  (=37.47  MHz) 


For  n  ■  6,  the  poles  in  radians/sec  are 
(-0.005  +  j 0.083)  x  109 
(-0.034  +  jO. 221)  x  109 
(-0.071  +  JO. 406)  x  109 


(=13.23  MHz) 
(=35.59  MHz) 
(=  65.9  MHz) 
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Here  a  good  approximation  to  the  poles  has  been  obtained  with  only 
four  poles.  Also,  there  seems  to  be  a  good  agreement  in  the  pole 
locations  obtained  from  the  original  integrated  data  and  the  numerically 
differentiated  data.  It  is  also  interesting  to  observe  that  indeed 
the  poles  are  occurring  approximately  at  odd  harmonics  of  the  fundamental. 

In  figure  2,  the  true  numerically  differentiated  data  are  plotted 
against  the  reconstructed  response  of  a  sixth  order  system.  The 
plot  has  been  normalized  to  unity  amplitude.  It  is  interesting  to 
note  that  there  is  a  close  agreement  even  in  the  very  early  times  of 
the  two  waveforms. 

8.  Conclusion 

The  pencil  of  function  method  has  been  applied  to  the  approximation 
of  an  arbitrary  function  by  a  sum  of  complex  exponentials.  The 
important  features  of  this  technique  are  the  natural  insensitivity 
to  noise  and  the  continuous  dependence  of  the  pole  locations  to  the 
integrated  square  error.  Finally  examples  have  been  presented  to 
illustrate  the  stability  of  the  pole  locations  yielded  by  the  pencil 
of  function  method. 
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10.  Appendix 

Let  A  and  li  be  two  square  NXN  matrices.  Let  the  columns  of 

A  and  B  be  represented  by  A^  and  for  i=l,2,...,N.  Denote 

C1  »  A,  if  v  =1 
Vi  1  1 

=  B±  if  v^=0 

and  the  matrix  constituted  by  the  columns  of 

C1  ,  C2  ,  C3  ,  ....  CN  as  C 
V1  V2  V3  VN  W"‘,VN 

Then 


det | [ A+B] |  =  l  det|c 
m=0 


-a>1,v2,..  .  ,vN 


where  in  the  second  summation  exactly  m  of  the  v^'s  equal  1  and  the 


rest  equal  0. 

This  is  a  standard  result  of  determinant  expansion. 


J 


Figure  2.  True  Response  Vs.  Reconstructed  Response  of  a  Sixth  Order  System  for  a  10m  long 
lm  Diameter  Conducting  Pipe. 


